#' Bounded Least Squares and Bounded Squared Hinge Loss Fcuntion Based
#' Twin Support Vector Machine
#'
#' \code{hinge_tsvm} is an R implementation of Hinge-TSVM
#'
#' @author Zhang Jiaqi.
#' @param X,y dataset and label.
#' @param C1,C2 plenty term.
#' @param kernel kernel function. The definitions of various kernel functions are as follows:
#' \describe{
#'     \item{linear:}{\eqn{u'v}{u'*v}}
#'     \item{poly:}{\eqn{(\gamma u'v + coef0)^{degree}}{(gamma*u'*v + coef0)^degree}}
#'     \item{rbf:}{\eqn{e^{(-\gamma |u-v|^2)}}{exp(-gamma*|u-v|^2)}}
#' }
#' @param gamma parameter for \code{'rbf'} and \code{'poly'} kernel. Default \code{gamma = 1/ncol(X)}.
#' @param degree parameter for polynomial kernel, default: \code{degree = 3}.
#' @param coef0 parameter for polynomial kernel,  default: \code{coef0 = 0}.
#' @param eps the precision of the optimization algorithm.
#' @param max.steps the number of iterations to solve the optimization problem.
#' @param solver \code{"dual"} is available.
#' @param fit_intercept if set \code{fit_intercept = TRUE},
#'                      the function will evaluates intercept.
#' @param randx parameter for reduce SVM, default \code{randx = 0.1}.
#' @return return \code{TSVMClassifier} object.
#' @export
bls_bsh_tsvm <- function(X, y, C1 = 1, C2 = C1, 
                         kernel = c("linear", "rbf", "poly"),
                         gamma = 1 / ncol(X), degree = 3, coef0 = 0,
                         lambda1 = 1, lambda2 = lambda1,
                         eps = 1e-5, eps.cccp = 1e-2,
                         max.steps = 4000, cccp.steps = 10,
                         solver = c("dual"), fit_intercept = TRUE,
                         randx = 1) {
  bsh_bls_tsvm_dual_solver <- function(KernelX, idx, C1, C2, lambda1, lambda2,
                                       eps, eps.cccp, 
                                       max.steps, cccp.steps) {
    update_delta_k <- function(z, lambda) {
      delta_k <- matrix(0, length(z))
      delta_k <- -z + z / ((1 + lambda*z^2)^2)
      return(as.matrix(delta_k))
    }
    update_sigma_k <- function(z, lambda) {
      sigma_k <- matrix(0, length(z))
      idx <- which(z >= 0)
      sigma_k[idx] <- -z[idx] + z[idx] / ((1 + lambda*z[idx]^2)^2)
      return(as.matrix(sigma_k))
    }
    
    X1 <- KernelX[idx, ]
    X2 <- KernelX[-idx, ]
    xn <- nrow(KernelX)
    xp <- ncol(KernelX)
    X1n <- nrow(X1)
    X2n <- xn - X1n
    
    solveres1 <- cholsolve(t(X1) %*% X1 + diag(1e-7, xp), t(KernelX))
    invX1TX1_X1T <- solveres1[, idx]
    invX1TX1_X2T <- solveres1[, -idx]
    X2_invX1TX1_X2T <- X2 %*% invX1TX1_X2T
    X1_invX1TX1_X2T <- X1 %*% invX1TX1_X2T
    
    solveres2 <- cholsolve(t(X2) %*% X2 + diag(1e-7, xp), t(KernelX))
    invX2TX2_X1T <- solveres2[, idx]
    invX2TX2_X2T <- solveres2[, -idx]
    X1_invX2TX2_X1T <- X1 %*% invX2TX2_X1T
    X2_invX2TX2_X1T <- X2 %*% invX2TX2_X1T
    
    dualH1  <- (X2_invX1TX1_X2T + diag(1/(C1), X2n))
    e2      <- matrix(1, X2n)
    duallb1 <- matrix(0, X2n)
    dualub1 <- matrix(Inf, X2n)
    
    dualH2 <- (X1_invX2TX2_X1T + diag(1/(C2), X1n))
    e1 <- matrix(1, X1n)
    duallb2 <- matrix(0, X1n)
    dualub2 <- matrix(Inf, X1n)
    
    a01 <- matrix(0, X2n)
    a02 <- matrix(0, X1n)
    
    delta_k1 <- matrix(0, X1n, 1)
    delta_k2 <- matrix(0, X2n, 1)
    
    sigma_k1 <- matrix(0, X2n, 1)
    sigma_k2 <- matrix(0, X1n, 1)
    
    u01 <- matrix(0, xp)
    u02 <- matrix(0, xp)
    
    for (i in 1:cccp.steps) {
      dualq1 <- -X2_invX1TX1_X2T %*% (sigma_k1*C1) -
                t(X1_invX1TX1_X2T) %*% (delta_k1) + e2
      a1 <- SupportVectorLab:::clip_dcd_optimizer(dualH1, dualq1, duallb1, dualub1,
                               eps, max.steps, a01)$x
      u1 <- -(invX1TX1_X2T %*% (a1 + C1*sigma_k1) + invX1TX1_X1T %*% delta_k1)
      f <- KernelX %*% u1
      delta_k1 <- update_delta_k(f[idx], lambda1)
      sigma_k1 <- update_sigma_k(1 + f[-idx], lambda1)
      if (norm(u1 - u01, type = "2") < eps.cccp) {
        break
      } else {
        u01 <- u1
        a01 <- a1
      }
    }
    
    for (i in 1:cccp.steps) {
      dualq2 <- -X1_invX2TX2_X1T %*% (sigma_k2*C2) +
                t(X2_invX2TX2_X1T) %*% delta_k2 + e1
      a2 <- SupportVectorLab:::clip_dcd_optimizer(dualH2, dualq2, duallb2, dualub2,
                               eps, max.steps, a02)$x
      u2 <- (invX2TX2_X1T %*% (a2 + C2*sigma_k2) - invX2TX2_X2T %*% delta_k2)
      f <- KernelX %*% u2
      delta_k2 <- update_delta_k(f[-idx], lambda2)
      sigma_k2 <- update_sigma_k(1 - f[idx], lambda2)
      if (norm(u2 - u02, type = "2") < eps.cccp) {
        break
      } else {
        u02 <- u2
        a02 <- a2
      }
    }
    BaseDualBLSBSHTSVMClassifier <- list("coef1" = as.matrix(u1),
                                         "coef2" = as.matrix(u2))
    return(BaseDualBLSBSHTSVMClassifier)
  }
  X <- as.matrix(X)
  y <- as.matrix(y)
  class_set <- sort(unique(y))
  idx <- which(y == class_set[1])
  y[idx] <- -1
  y[-idx] <- 1
  y <- as.matrix(as.numeric(y))
  if (length(class_set) > 2) {
    stop("The number of class should less 2!")
  }
  kernel <- match.arg(kernel)
  solver <- match.arg(solver)
  kso <- SupportVectorLab:::kernel_select_option(X, kernel, "primal", randx,
                                         gamma, degree, coef0)
  KernelX <- kso$KernelX
  Kw <- kso$KernelX[kso$sample_idx, ]
  X <- kso$X
  if (fit_intercept == TRUE) {
    KernelX <- cbind(KernelX, 1)
  }
  if (solver == "dual") {
    solver.res <- bsh_bls_tsvm_dual_solver(KernelX, idx, C1, C2,
                                           lambda1, lambda2,
                                           eps, eps.cccp,
                                           max.steps, cccp.steps)
  }
  TSVMClassifier <- list("X" = X, "y" = y, "class_set" = class_set,
                         "C1" = C1, "C2" = C2, "kernel" = kernel,
                         "gamma" = gamma, "degree" = degree, "coef0" = coef0,
                         "solver" = solver, "coef1" = solver.res$coef1,
                         "coef2" = solver.res$coef2,
                         "fit_intercept" = fit_intercept,
                         "Kw" = Kw)
  class(TSVMClassifier) <- "TSVMClassifier"
  return(TSVMClassifier)
}